options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)

execute mcmc sampling

goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
  mcmc=mdl$sample(
  data=data,
  seed=1,
  chains=4,
  iter_sampling=smp,
  iter_warmup=wrm,
  thin=th,
  refresh=1000
  )
  mcmc
}

see mcmc result and parameters

seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
  print(mcmc)
  prs=mcmc$metadata()$model_params[-1] # reject lp__
  for(pr in prs){
    if(grepl('^y',pr)) next # not show predictive value "y*" information
    if(exc!='' && grepl(paste0('^',exc),pr)) next
    drw=mcmc$draws(pr)
    if(ch){
      par(mfrow=c(4,1),mar=c(1,5,1,1))
      drw[,1,] |> plot(type='l',xlab='',ylab=pr)
      drw[,2,] |> plot(type='l',xlab='',ylab=pr)
      drw[,3,] |> plot(type='l',xlab='',ylab=pr)
      drw[,4,] |> plot(type='l',xlab='',ylab=pr)
      par(mar=c(3,5,3,3))
    }

    par(mfrow=c(1,2))
    drw |> hist(main=pr,xlab='')
    drw |> density() |> plot(main=pr)    
  }
}

stan basic

normal distribution

ex1-0.stan

normal distribution

data{
  int N;
  vector[N] y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~normal(m,s);
}
generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
mdl=cmdstan_model('./ex1-0.stan')
y=rnorm(10,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.66    1.74    2.20    2.17    2.40    3.42
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')  

data=list(N=length(y),y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -32.4689 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        9      -2.00847   0.000241842   0.000187564           1           1       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.3 seconds.
mle
##  variable estimate
##      lp__    -2.01
##      m        2.17
##      s        0.74
##      y1       1.22
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
mcmc
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -3.35  -3.01 1.04 0.75 -5.46 -2.34 1.00      803      925
##      m     2.16   2.16 0.30 0.28  1.69  2.62 1.00     1091      918
##      s     0.91   0.86 0.25 0.21  0.60  1.37 1.00     1380     1165
##      y1    2.14   2.15 0.99 0.91  0.55  3.72 1.00     1907     1879
mcmc$metadata()$stan_variables
## [1] "lp__" "m"    "s"    "y1"
mcmc$metadata()$model_params
## [1] "lp__" "m"    "s"    "y1"
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -3.35  -3.01 1.04 0.75 -5.46 -2.34 1.00      803      925
##      m     2.16   2.16 0.30 0.28  1.69  2.62 1.00     1091      918
##      s     0.91   0.86 0.25 0.21  0.60  1.37 1.00     1380     1165
##      y1    2.14   2.15 0.99 0.91  0.55  3.72 1.00     1907     1879

poisson distribution

y=rpois(20,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.75    2.50    2.95    3.25   10.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=length(y),y=y)

ex2-1.stan

in case fit poisson dist. to normal dist.

data{
  int N;
  vector[N] y;
}

parameters{
  real<lower=0> m;
  real<lower=0> s;
}

model{
  y~normal(m,s);
}

generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
# when fitting poisson dist. sample to normal dist.
mdl=cmdstan_model('./ex2-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -33.8038 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -26.7663   0.000851482    0.00122998           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -26.77
##      m        2.95
##      s        2.31
##      y1       8.05
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -25.88 -25.59 1.02 0.77 -27.97 -24.88 1.00      786      807
##      m      2.96   2.94 0.58 0.58   2.04   3.91 1.00     1065      870
##      s      2.53   2.48 0.43 0.40   1.93   3.31 1.00     1610     1257
##      y1     2.96   2.88 2.63 2.52  -1.43   7.37 1.00     1932     1974

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1')     

ex2-2.stan

fit to poisson dist.

data{
  int N;
  array[N] int y;
}

parameters{
  real<lower=0> l;
}

model{
  y~poisson(l);
}

generated quantities{
  int y1;
  y1=poisson_rng(l);
}
mdl=cmdstan_model('./ex2-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -20.7726 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       4.82651    0.00029093   0.000121597           1           1        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__     4.83
##      l        2.95
##      y1       4.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##      lp__ 5.44   5.71 0.66 0.28 4.05 5.92 1.00      833     1284
##      l    3.02   3.00 0.38 0.37 2.39 3.67 1.00      819      938
##      y1   3.00   3.00 1.79 1.48 0.00 6.00 1.00     1710     1869

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

difference test

ex3.stan

mean difference of 2 normal distributions

data{
  int N;
  vector[N] a;
  vector[N] b;
}
parameters{
  real ma;
  real<lower=0> sa;
  real mb;
  real<lower=0> sb;
}
model{
  a~normal(ma,sa);
  b~normal(mb,sb);
}
generated quantities{
  real d;
  d=mb-ma;
}
n=20
a=rnorm(n,10,1)
b=rnorm(n,11,2)
data=list(N=n,a=a,b=b)

mdl=cmdstan_model('./ex3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -33482.6 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       24      -33.3995   0.000938568    0.00156273      0.9537      0.9537       59    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -33.40
##      ma      10.14
##      sa       0.78
##      mb      11.13
##      sb       2.49
##      d        0.99
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -34.75 -34.42 1.42 1.26 -37.52 -33.03 1.00      799     1160
##      ma    10.15  10.15 0.19 0.20   9.83  10.45 1.00     1892     1204
##      sa     0.87   0.84 0.16 0.14   0.65   1.17 1.00     2124     1341
##      mb    11.12  11.13 0.59 0.56  10.10  12.08 1.01      476      619
##      sb     2.74   2.67 0.49 0.46   2.08   3.67 1.00     2226     1544
##      d      0.97   0.98 0.62 0.60  -0.06   1.98 1.00      537      709

d=mcmc$draws('d')
d
## # A draws_array: 500 iterations, 4 chains, and 1 variables
## , , variable = d
## 
##          chain
## iteration     1    2     3    4
##         1 2.215 0.92  1.68 0.76
##         2 0.068 0.55  2.43 0.98
##         3 0.659 0.80  1.63 1.19
##         4 1.960 0.47  1.51 2.31
##         5 1.953 0.39 -0.41 2.17
## 
## # ... with 495 more iterations
mean(d>0)
## [1] 0.935

normal regression

ex4-1.stan

single normal regression

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}
# make prediction with R
n=20
x=runif(n,0,100)
y=rnorm(n,x*3+10,10)
par(mfrow=c(1,1))
plot(x,y)

data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -50125.8 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
## Error evaluating model log probability: Non-finite gradient. 
## Error evaluating model log probability: Non-finite gradient. 
##       34      -58.8763    0.00241247   0.000958411       0.945       0.945       90    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -58.88
##      b0      15.00
##      b1       2.95
##      s       11.52
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -58.16 -57.66 1.55 1.09 -61.15 -56.62 1.01      302      225
##      b0    15.18  15.04 6.50 5.60   5.25  26.29 1.02      176      138
##      b1     2.95   2.95 0.10 0.09   2.78   3.11 1.01      229      241
##      s     13.22  12.81 2.47 2.16   9.97  17.75 1.01      746      594

b0=mcmc$draws('b0')
b1=mcmc$draws('b1')
b=tibble(b0=c(b0),b1=c(b1)) |> sample_n(100)
x1=runif(nrow(b),min(x),max(x))
xy=tibble(x=x1,y=b$b0+b$b1*x1)
par(mfrow=c(1,1))
plot(xy)

ex4-2.stan

single normal regression, prediction from new explanatory variable x1

data{
  int N;
  vector[N] x;
  vector[N] y;
  int N1;
  vector[N1] x1;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N1] m;
  vector[N1] y1;
  for(i in 1:N1){
    m[i]=b0+b1*x1[i];
    y1[i]=normal_rng(m[i],s);
  }
}
# make prediction with stan
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition

data=list(N=n,x=x,y=y,N1=n1,x1=x1)

mdl=cmdstan_model('./ex4-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -27844.4 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       27      -58.8763     0.0076696    0.00260617           1           1       71    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##    lp__     -58.88
##    b0        14.99
##    b1         2.95
##    s         11.52
##    m[1]      22.35
##    m[2]      52.97
##    m[3]      83.60
##    m[4]     114.23
##    m[5]     144.85
##    m[6]     175.48
##    m[7]     206.11
##    m[8]     236.74
##    m[9]     267.36
##    m[10]    297.99
##    y1[1]     13.58
##    y1[2]     49.62
##    y1[3]     88.66
##    y1[4]    124.27
##    y1[5]    151.62
##    y1[6]    157.93
##    y1[7]    206.78
##    y1[8]    227.16
##    y1[9]    303.72
##    y1[10]   303.28
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
mcmc$metadata()$stan_variables
## [1] "lp__" "b0"   "b1"   "s"    "m"    "y1"
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -58.04 -57.76  1.29  1.09 -60.57 -56.62 1.01      445      595
##    b0      14.69  14.65  5.61  5.40   5.48  23.77 1.02      249      180
##    b1       2.95   2.95  0.09  0.09   2.80   3.10 1.01      326      337
##    s       13.07  12.79  2.40  2.21   9.93  17.33 1.00     1100      881
##    m[1]    22.05  22.04  5.41  5.22  13.20  30.83 1.02      250      183
##    m[2]    52.72  52.68  4.64  4.38  45.01  60.25 1.02      253      185
##    m[3]    83.40  83.34  3.95  3.80  76.81  89.93 1.02      268      210
##    m[4]   114.07 113.99  3.39  3.28 108.48 119.63 1.01      313      323
##    m[5]   144.74 144.68  3.03  2.90 139.85 149.74 1.00      450      580
##    m[6]   175.42 175.39  2.95  2.84 170.64 180.24 1.00      967      924
##    m[7]   206.09 206.07  3.17  3.08 200.69 211.20 1.00     2319     1649
##    m[8]   236.76 236.72  3.64  3.38 230.58 242.73 1.01     2458     1416
##    m[9]   267.43 267.36  4.27  3.95 260.30 274.49 1.01     1869     1405
##    m[10]  298.11 298.05  5.01  4.82 289.95 306.33 1.00     1221     1348
##    y1[1]   21.81  21.85 14.24 13.78  -1.53  44.47 1.00      979     1344
##    y1[2]   52.40  52.33 14.02 13.50  29.62  75.10 1.00     1010     1136
##    y1[3]   83.15  83.12 13.82 13.27  60.58 104.87 1.00     1270     1614
##    y1[4]  114.63 114.20 13.73 13.61  93.29 137.27 1.00     1617     1858
##    y1[5]  144.31 144.88 13.76 13.42 121.85 166.26 1.00     1833     1720
##    y1[6]  175.56 175.58 13.52 12.70 153.84 197.48 1.00     1873     2031
##    y1[7]  206.43 206.31 13.81 12.91 183.85 228.30 1.00     1800     1984
##    y1[8]  236.46 236.72 13.80 13.16 213.56 259.08 1.00     1823     1855
##    y1[9]  266.85 266.52 13.97 13.34 244.01 290.08 1.00     1823     1741
##    y1[10] 298.22 298.30 14.38 13.73 274.54 321.91 1.00     2153     1598

m=mcmc$draws('m')
summary(m)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 m[1]      22.1   22.0  5.41  5.22  13.2  30.8  1.02     250.     184.
##  2 m[2]      52.7   52.7  4.64  4.38  45.0  60.3  1.02     254.     186.
##  3 m[3]      83.4   83.3  3.95  3.80  76.8  89.9  1.02     269.     211.
##  4 m[4]     114.   114.   3.39  3.28 108.  120.   1.01     313.     324.
##  5 m[5]     145.   145.   3.03  2.90 140.  150.   1.00     451.     580.
##  6 m[6]     175.   175.   2.95  2.84 171.  180.   1.00     968.     924.
##  7 m[7]     206.   206.   3.17  3.08 201.  211.   1.00    2319.    1649.
##  8 m[8]     237.   237.   3.64  3.38 231.  243.   1.01    2459.    1417.
##  9 m[9]     267.   267.   4.27  3.95 260.  274.   1.01    1870.    1406.
## 10 m[10]    298.   298.   5.01  4.82 290.  306.   1.00    1221.    1349.
smm=summary(m)

y1=mcmc$draws('y1')
summary(y1)
## # A tibble: 10 × 10
##    variable  mean median    sd   mad     q5   q95  rhat ess_bulk ess_tail
##    <chr>    <dbl>  <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
##  1 y1[1]     21.8   21.8  14.2  13.8  -1.53  44.5 1.00      979.    1344.
##  2 y1[2]     52.4   52.3  14.0  13.5  29.6   75.1 1.00     1011.    1137.
##  3 y1[3]     83.2   83.1  13.8  13.3  60.6  105.  1.00     1270.    1615.
##  4 y1[4]    115.   114.   13.7  13.6  93.3  137.  0.999    1617.    1859.
##  5 y1[5]    144.   145.   13.8  13.4 122.   166.  1.00     1833.    1720.
##  6 y1[6]    176.   176.   13.5  12.7 154.   197.  1.00     1873.    2032.
##  7 y1[7]    206.   206.   13.8  12.9 184.   228.  1.00     1800.    1985.
##  8 y1[8]    236.   237.   13.8  13.2 214.   259.  1.00     1823.    1855.
##  9 y1[9]    267.   267.   14.0  13.3 244.   290.  1.00     1823.    1741.
## 10 y1[10]   298.   298.   14.4  13.7 275.   322.  1.00     2153.    1599.
smy=summary(y1)

xy=tibble(x=x1,m=smm$median,ml5=smm$q5,mu5=smm$q95,yl5=smy$q5,yu5=smy$q95)

par(mfrow=c(1,1))
xlim=c(min(x),max(x))
ylim=c(min(y),max(y))
plot(x,y,
     xlim=xlim,ylim=ylim, xlab='x',ylab='y',col='red')
par(new=T)
plot(xy$x,xy$m,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='black')
par(new=T)
plot(xy$x,xy$ml5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$mu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='darkgray')
par(new=T)
plot(xy$x,xy$yl5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')
par(new=T)
plot(xy$x,xy$yu5,type='l',
     xlim=xlim,ylim=ylim, xlab='',ylab='',col='lightgray')

qplot(x,y,col=I('red'))+
  geom_line(aes(x=x,y=m),xy,col='black')+
  geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
  geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
  geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
  geom_line(aes(x=x,y=yu5),xy,col='lightgray')

ex4-3.stan

single normal regression, prediction from original explanatory variable x

data{
  int N;
  vector[N] x;
  vector[N] y;
}

parameters{
  real b0;
  real b1;
  real<lower=0> s;
}

model{
  y~normal(b0+b1*x,s);
}

generated quantities{
  vector[N] m;
  vector[N] y1;
  for(i in 1:N){
    m[i]=b0+b1*x[i];
    y1[i]=normal_rng(m[i],s);
  }
}
data=list(N=n,x=x,y=y)

mdl=cmdstan_model('./ex4-3.stan')

mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc,'m',ch=F)
##  variable   mean median    sd   mad     q5    q95 rhat ess_bulk ess_tail
##    lp__   -58.15 -57.84  1.38  1.19 -60.76 -56.62 1.01      281      272
##    b0      14.44  15.07  6.47  6.43   3.59  24.17 1.01      167      143
##    b1       2.96   2.95  0.10  0.10   2.80   3.13 1.01      219      216
##    s       13.13  12.83  2.40  2.27   9.75  17.57 1.01      808      841
##    m[1]   146.73 146.95  3.18  3.11 141.35 151.75 1.01      279      258
##    m[2]   265.58 265.60  4.12  3.90 258.85 272.30 1.00     1452     1279
##    m[3]   200.49 200.59  3.01  2.75 195.28 205.48 1.00     1535      901
##    m[4]   215.81 215.85  3.16  2.93 210.44 221.07 1.00     2344     1175
##    m[5]    26.49  27.10  6.10  6.07  16.33  35.67 1.01      167      136
##    m[6]   278.42 278.53  4.44  4.21 271.27 285.67 1.00     1169     1150
##    m[7]   298.14 298.27  4.97  4.78 289.99 306.24 1.00      907     1082
##    m[8]   151.22 151.43  3.12  3.03 145.95 156.18 1.01      301      275
##    m[9]    44.95  45.51  5.55  5.54  35.59  53.21 1.01      167      133
##    m[10]  150.42 150.63  3.13  3.04 145.12 155.40 1.01      296      275
##    m[11]  172.14 172.24  2.96  2.84 167.12 176.99 1.00      511      507
##    m[12]  272.57 272.62  4.29  4.07 265.63 279.59 1.00     1273     1212
##    m[13]  101.53 101.88  4.04  3.96  94.47 107.61 1.01      186      160
##    m[14]  131.78 132.06  3.41  3.28 125.96 137.03 1.01      229      223
##    m[15]  281.29 281.42  4.51  4.29 273.98 288.69 1.00     1123     1065
##    m[16]   84.90  85.34  4.45  4.35  77.24  91.54 1.01      177      144
##    m[17]  282.64 282.78  4.55  4.34 275.27 290.10 1.00     1102     1065
##    m[18]  296.20 296.35  4.91  4.72 288.15 304.24 1.00      928     1073
##    m[19]   21.81  22.44  6.24  6.21  11.40  31.25 1.01      167      139
##    m[20]   24.46  25.09  6.16  6.11  14.17  33.76 1.01      167      138
##    y1[1]  146.81 147.07 13.85 13.23 123.78 169.41 1.01     1842     1850
##    y1[2]  265.51 265.72 13.71 13.29 242.70 288.06 1.00     1979     1794
##    y1[3]  200.03 200.42 13.99 13.05 177.31 222.05 1.00     1879     1825
##    y1[4]  216.19 216.34 13.97 13.53 193.34 238.59 1.00     1912     1848
##    y1[5]   26.84  27.16 14.87 14.42   2.60  50.41 1.00     1094     1540
##    y1[6]  279.31 278.90 13.97 13.30 257.20 302.44 1.00     1848     1988
##    y1[7]  298.62 298.46 14.27 13.26 275.75 321.73 1.00     2039     1635
##    y1[8]  151.51 151.78 13.79 13.20 128.37 173.46 1.00     1876     1761
##    y1[9]   45.56  45.56 14.57 14.60  21.29  68.85 1.00      751     1492
##    y1[10] 150.44 150.36 13.22 12.44 128.80 172.38 1.00     2018     1773
##    y1[11] 171.94 172.02 13.37 12.94 149.26 193.84 1.00     2187     1957
##    y1[12] 272.22 271.98 14.35 13.70 248.83 295.12 1.00     1955     1895
##    y1[13] 101.34 101.77 13.72 13.06  78.58 122.87 1.00     1590     1813
##    y1[14] 132.00 132.33 13.88 13.23 108.43 154.81 1.00     2008     1820
##    y1[15] 281.55 281.64 13.78 12.87 258.73 303.51 1.00     2025     1673
##    y1[16]  84.90  85.07 14.42 13.72  61.50 107.17 1.00     1192     1370
##    y1[17] 283.34 283.04 14.07 12.95 260.59 307.26 1.00     1861     1897
##    y1[18] 295.94 295.57 14.26 14.04 272.46 319.17 1.00     1768     1933
##    y1[19]  21.79  21.80 14.73 14.31  -3.10  45.41 1.00      863     1135
##    y1[20]  24.35  24.63 14.85 14.23   0.00  47.42 1.00      645      895

y1=mcmc$draws('y1')
smy=summary(y1)

par(mfrow=c(1,1))
plot(y,smy$median,xlab='obs.',ylab='prd.')
abline(0,1)

qplot(y,smy$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)

par(mfrow=c(1,2))
hist(y-smy$median,xlab='obs.-prd.',main='residual')
density(y-smy$median) |> plot(xlab='obs.-prd.',main='residual')

estimate correlation

ex4-41.stan

use covariance matrix and multinormal distribution

data{
  int N;
  array[N] vector[2] y;
}

parameters{
  vector[2] m;
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  y~multi_normal(m,cv);
}

ex4-42.stan

use covariance matrix and wishart distribution

data{
  int N;
  matrix[2,2] dp;
}

parameters{
  vector<lower=0>[2] s;
  real<lower=-1,upper=1> r;
}

transformed parameters{
  cov_matrix[2] cv;
  cv[1,1]=s[1]^2;
  cv[2,2]=s[2]^2;
  cv[1,2]=r*s[1]*s[2];
  cv[2,1]=r*s[1]*s[2];
}

model{
  dp~wishart(N-1,cv);
}

ex4-43.stan

use correlation matrix and wishart distribution

data{
  int N;
  int M;
  matrix[M,M] dp;
}

parameters{
  vector<lower=0>[M] s;
  corr_matrix[M] r;
}

transformed parameters{
  cov_matrix[M] cv;
  cv=quad_form_diag(r,s);
}

model{
  dp~wishart(N-1,cv);
}

ex4-44.stan use covariance matrix and multinormal distribution

data{
  int N;
  int K;
  array[N] vector[K] y;
}
parameters{
  vector[K] m;
  cov_matrix[K] cov;
}
model{
  y~multi_normal(m,cov);
}
library(MASS)
n=20
y=mvrnorm(n,c(0,0),matrix(c(1,0.7,0.7,1),2),2)
cov(y)
##       [,1]  [,2]
## [1,] 1.213 0.674
## [2,] 0.674 0.902
cor(y)
##       [,1]  [,2]
## [1,] 1.000 0.644
## [2,] 0.644 1.000
plot(y)

#estimate covariance
mdl=cmdstan_model('./ex4-41.stan')

data=list(N=n,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -132.14 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       23      -14.5172   5.04246e-05    0.00107316      0.3566           1       28    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -14.52
##   m[1]       -0.02
##   m[2]        0.24
##   s[1]        1.07
##   s[2]        0.93
##   r           0.64
##   cv[1,1]     1.15
##   cv[2,1]     0.64
##   cv[1,2]     0.64
##   cv[2,2]     0.86
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -18.41 -18.02 1.73 1.53 -21.71 -16.35 1.00      812     1202
##   m[1]     -0.03  -0.04 0.26 0.25  -0.45   0.43 1.01     1054     1071
##   m[2]      0.24   0.24 0.24 0.22  -0.14   0.62 1.00     1276     1050
##   s[1]      1.18   1.16 0.21 0.19   0.89   1.56 1.00     1305     1236
##   s[2]      1.03   1.00 0.19 0.17   0.78   1.36 1.00     1333     1505
##   r         0.58   0.61 0.16 0.15   0.26   0.80 1.01      427      660
##   cv[1,1]   1.44   1.34 0.55 0.44   0.80   2.45 1.00     1305     1236
##   cv[2,1]   0.75   0.68 0.39 0.32   0.25   1.42 1.01      517      791
##   cv[1,2]   0.75   0.68 0.39 0.32   0.25   1.42 1.01      517      791
##   cv[2,2]   1.10   1.01 0.44 0.35   0.61   1.85 1.00     1333     1505

#estimate correlation,covariance
#deviation product sum matirix folows Wishart dist.
#deviation product sum = covariance * n-1 or n

mdl=cmdstan_model('./ex4-42.stan')

data=list(N=n,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -26.1847 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -14.7659   0.000218707   1.49717e-05           1           1       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -14.77
##   s[1]        1.10
##   s[2]        0.95
##   r           0.64
##   cv[1,1]     1.21
##   cv[2,1]     0.67
##   cv[1,2]     0.67
##   cv[2,2]     0.90
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -17.44 -17.14 1.22 1.03 -19.81 -16.11 1.00      755     1168
##   s[1]      1.19   1.16 0.21 0.19   0.91   1.59 1.00     1223     1318
##   s[2]      1.03   1.00 0.19 0.17   0.77   1.37 1.00     1221     1190
##   r         0.60   0.61 0.15 0.14   0.33   0.82 1.01      532      721
##   cv[1,1]   1.47   1.35 0.56 0.44   0.83   2.52 1.00     1223     1318
##   cv[2,1]   0.77   0.69 0.39 0.31   0.30   1.53 1.01      596      725
##   cv[1,2]   0.77   0.69 0.39 0.31   0.30   1.53 1.01      596      725
##   cv[2,2]   1.09   1.00 0.42 0.35   0.59   1.89 1.00     1221     1190

#estimate correlation,covariance
mdl=cmdstan_model('./ex4-43.stan')

m=2
data=list(N=n,M=m,dp=cov(y)*(n-1))

mle=mdl$optimize(data=data)
## Initial log joint probability = -99.3293 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -14.7659   0.000171539   0.000536133           1           1       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##   lp__      -14.77
##   s[1]        1.10
##   s[2]        0.95
##   r[1,1]      1.00
##   r[2,1]      0.64
##   r[1,2]      0.64
##   r[2,2]      1.00
##   cv[1,1]     1.21
##   cv[2,1]     0.67
##   cv[1,2]     0.67
##   cv[2,2]     0.90
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##   lp__    -16.70 -16.46 1.15 1.03 -18.90 -15.39 1.01      921     1403
##   s[1]      1.18   1.16 0.20 0.19   0.90   1.56 1.00     1256     1360
##   s[2]      1.02   0.99 0.17 0.16   0.78   1.32 1.01     1189     1214
##   r[1,1]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   r[2,1]    0.59   0.61 0.15 0.15   0.31   0.81 1.00      889     1021
##   r[1,2]    0.59   0.61 0.15 0.15   0.31   0.81 1.00      889     1021
##   r[2,2]    1.00   1.00 0.00 0.00   1.00   1.00   NA       NA       NA
##   cv[1,1]   1.44   1.33 0.52 0.44   0.80   2.44 1.00     1256     1360
##   cv[2,1]   0.74   0.68 0.36 0.30   0.29   1.41 1.00      831      764
##   cv[1,2]   0.74   0.68 0.36 0.30   0.29   1.41 1.00      831      764
##   cv[2,2]   1.06   0.98 0.38 0.31   0.60   1.76 1.01     1189     1214

mdl=cmdstan_model('./ex4-44.stan')

k=2
data=list(N=n,K=k,y=y)

mle=mdl$optimize(data=data)
## Initial log joint probability = -1569.73 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       16      -14.5172   8.20179e-05   0.000961096      0.8987      0.8987       21    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##  lp__       -14.52
##  m[1]        -0.02
##  m[2]         0.24
##  cov[1,1]     1.15
##  cov[2,1]     0.64
##  cov[1,2]     0.64
##  cov[2,2]     0.86
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc,ch=F)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##  lp__     -16.13 -15.79 1.67 1.49 -19.28 -14.09 1.00      815     1069
##  m[1]      -0.03  -0.02 0.29 0.27  -0.53   0.43 1.01      846      817
##  m[2]       0.24   0.23 0.25 0.25  -0.18   0.65 1.01     1132      942
##  cov[1,1]   1.78   1.62 0.78 0.60   0.92   3.14 1.00     1522     1339
##  cov[2,1]   0.98   0.88 0.51 0.41   0.39   1.93 1.00     1460     1146
##  cov[1,2]   0.98   0.88 0.51 0.41   0.39   1.93 1.00     1460     1146
##  cov[2,2]   1.29   1.17 0.53 0.41   0.69   2.29 1.00     1406     1191

distributions

normal dist.

distribution of mean, it's from central limit theorem
y~N(m,s), y(-Inf,Inf), E[y]=m, V[y]=s^2

ex1-0.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real m;
  real<lower=0> s;
}
model{
  y~normal(m,s);
}
generated quantities{
  real y1;
  y1=normal_rng(m,s);
}
n=20
y=rnorm(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.029   0.835   1.536   1.473   2.105   2.867
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-0.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -38.7162 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       12      -4.99263   0.000103972   9.41093e-05      0.9847      0.9847       18    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -4.99
##      m        1.47
##      s        0.78
##      y1       2.39
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -6.22  -5.94 0.96 0.71 -8.18 -5.29 1.00      947     1127
##      m     1.48   1.48 0.18 0.18  1.18  1.78 1.00     1363     1058
##      s     0.84   0.83 0.14 0.14  0.64  1.10 1.00     1743     1208
##      y1    1.46   1.45 0.89 0.84  0.03  2.93 1.00     1856     1888

Bernoulli dist.

The event with probablity p occur y=1, not occur y=0 
y~Ber(p), y{0,1}, E[y]=p, V[y]=p(1-p)

ex1-1.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  y~bernoulli(p);
}
generated quantities{
  real s;
  s=(p*(1-p))^.5;
}
n=10
y=sample(0:1,n,replace=T)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     1.0     1.0     0.8     1.0     1.0
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-1.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -5.0084 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        3      -5.00402    0.00127944   3.77503e-05           1           1        6    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -5.00
##      p        0.80
##      s        0.40
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -7.27  -7.00 0.69 0.35 -8.76 -6.75 1.00      909      859
##      p     0.75   0.76 0.12 0.12  0.52  0.92 1.00      617      827
##      s     0.41   0.43 0.07 0.07  0.27  0.50 1.00      566      827

geometric dist.

The event with probablity p occur after y trials(failure) 
y~Ge(p), y[0,Inf), P(Y=y)=p(1−p)^y
E[y]=1/p, V[y]=(1-p)/p^2

ex1-2.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
}
model{
  for (n in 1:N) {
    target+=log(p)+y[n]*log(1-p);
  }
}
generated quantities{
  real m;
  m=1/p;
  real s;
  s=((1-p)/p^2)^.5;
}
n=20
y=rgeom(n,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    0.75    1.50    3.05    6.00   10.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-2.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -45.641 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        3      -45.2724    0.00360832   0.000134258      0.9482      0.9482        5    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -45.27
##      p        0.25
##      m        4.05
##      s        3.51
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -47.42 -47.17 0.65 0.30 -48.74 -46.95 1.00      903     1451
##      p      0.25   0.25 0.05 0.05   0.18   0.33 1.00      747      952
##      m      4.08   3.94 0.79 0.71   3.03   5.56 1.00      747      952
##      s      3.55   3.40 0.80 0.72   2.48   5.03 1.00      747      952

negative binomial dist.

The event with probablity p occur a times after y-a trials 
y~NB(a,p), y[0,Inf), P(Y=y)=p^a*(1−p)^(y-a)
E[y]=a*(1-p)/p, V[y]=a*(1-p)/p^2

ex1-3.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0,upper=1> p;
  real<lower=1> a;
}
model{
  y~neg_binomial(a,p);
}
generated quantities{
  real m;
  m=a*(1-p)/p;
  real s;
  s=(a*(1-p)/p^2)^.5;
  int y1;
  y1=neg_binomial_rng(a,p);
}
n=50
a=3
y=rnbinom(n,a,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    3.00    5.00    5.90    7.75   18.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -143.873 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11      -137.898    0.00382547   0.000802942      0.9311      0.9311       13    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -137.90
##      p        0.41
##      a        2.43
##      m        3.47
##      s        2.90
##      y1      13.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.3 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -139.93 -139.57 1.11 0.76 -142.12 -138.86 1.00      481      541
##      p       0.52    0.50 0.16 0.16    0.30    0.82 1.01      443      313
##      a       3.03    2.89 0.89 0.88    1.81    4.69 1.01      448      455
##      m       2.84    2.86 1.07 1.03    1.00    4.55 1.01      496      313
##      s       2.44    2.39 0.85 0.82    1.11    3.84 1.01      466      326
##      y1      5.78    5.00 4.34 4.45    1.00   14.00 1.00     1916     1799

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

binomial dist.

The event with probablity p occur y times on n trials
y~B(n,p), y{0,1...n}, E[y]=np, V[y]=np(1-p)

ex2-3.stan

data{
  int N;
  int n;
  array[N] int y;
}
parameters{
  real<lower=0,upper=0> p;
}
model{
  y~binomial(n,p);
}
generated quantities{
  real m;
  m=n*p;
  real s;
  s=(n*p*(1-p))^.5;
  int y1;
  y1=binomial_rng(n,p);
}
n=20
n0=10
y=rbinom(n,n0,0.3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     2.0     3.0     3.1     4.0     6.0
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex2-3.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -144.737 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5       -123.82    0.00185335   0.000182482           1           1        7    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__  -123.82
##      p        0.31
##      m        3.10
##      s        1.46
##      y1       2.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable    mean  median   sd  mad      q5     q95 rhat ess_bulk ess_tail
##      lp__ -125.86 -125.60 0.69 0.33 -127.22 -125.36 1.00     1168     1356
##      p       0.31    0.31 0.03 0.03    0.26    0.37 1.00      983     1175
##      m       3.14    3.13 0.33 0.34    2.61    3.69 1.00      983     1175
##      s       1.46    1.47 0.04 0.04    1.39    1.53 1.00      983     1175
##      y1      3.11    3.00 1.51 1.48    1.00    6.00 1.00     1747     1964

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

ex2-31.stan

trial n is varied from each sample

data{
  int N;
  array[N] int n;
  array[N] int y;
}
parameters{
  real<lower=0> p;
}
model{
  y~binomial(n,p);
}
n=20
n0=floor(runif(n,1,10))
y=rbinom(n,n0,0.3)

data=list(N=n,n=n0,y=y)

mdl=cmdstan_model('./ex2-31.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -81.2821 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        4       -66.326    0.00582635    0.00143967           1           1        6    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -66.33
##      p        0.29
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -68.37 -68.10 0.68 0.27 -69.76 -67.90 1.00     1104     1270
##      p      0.29   0.29 0.04 0.04   0.23   0.36 1.01      707      968

Poisson dist.

The event occur l times in unit range, the event occur y times. 
y~Po(l), y{0,1...Inf}, E[y]=l, V[y]=l

ex1-4.stan

data{
  int N;
  array[N] int y;
}
parameters{
  real<lower=0> l;
}
model{
  y~poisson(l);
}
generated quantities{
  real s;
  s=l^.5;
  int y1;
  y1=poisson_rng(l);
}
n=20
y=rpois(n,3)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    1.75    3.00    3.30    4.00   12.00
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-4.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -58.9394 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        5       12.7989    0.00515492   0.000519729      0.9927      0.9927        8    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    12.80
##      l        3.30
##      s        1.82
##      y1       1.00
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ 13.49  13.77 0.70 0.31 12.02 14.00 1.00     1172     1478
##      l     3.37   3.35 0.41 0.41  2.72  4.09 1.01      881     1161
##      s     1.83   1.83 0.11 0.11  1.65  2.02 1.01      881     1161
##      y1    3.36   3.00 1.87 1.48  1.00  7.00 1.00     1881     1892

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

exponential dist.

The event occur l times in unit range, the event take time y to occur
y~ex(l), y>0, E[y]=1/l, V[y]=1/l^2

ex1-5.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> l;
}
model{
  y~exponential(l);
}
generated quantities{
  real m;
  m=1/l;
  real s;
  s=1/l;
  real y1;
  y1=exponential_rng(l);
}
n=20
y=rexp(n,2)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.035   0.166   0.443   0.473   0.666   1.337
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-5.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -12.4359 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6      -5.00923   0.000141713   6.78712e-05      0.9789      0.9789        9    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__    -5.01
##      l        2.12
##      m        0.47
##      s        0.47
##      y1       2.30
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##      lp__ -4.74  -4.47 0.70 0.32 -6.21 -4.24 1.00     1026     1339
##      l     2.20   2.16 0.48 0.47  1.49  3.04 1.00      758     1014
##      m     0.48   0.46 0.11 0.10  0.33  0.67 1.00      758     1014
##      s     0.48   0.46 0.11 0.10  0.33  0.67 1.00      758     1014
##      y1    0.47   0.30 0.50 0.33  0.03  1.44 1.00     1755     1759

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

gamma dist.

The event occur l times in unit range, the event take time y to occur a times
y~Ga(a,l), y>0, E[y]=a/l, V[y]=a/l^2

ex1-6.stan

data{
  int N;
  vector[N] y;
}
parameters{
  real<lower=0> a;
  real<lower=0> l;
}
model{
  y~gamma(a,l);
}
n=20
y=rgamma(n,2,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.86    1.59    2.17    2.27    2.78    5.15
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-6.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -137.813 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       13      -26.5906   0.000411502   0.000416302        0.92        0.92       14    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -26.59
##      a        5.40
##      l        2.38
##      m        2.27
##      s        0.98
##      y1       1.42
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -24.84 -24.54 0.96 0.70 -26.84 -23.90 1.01      564      932
##      a      6.12   5.92 1.73 1.66   3.60   9.36 1.01      400      402
##      l      2.72   2.65 0.79 0.79   1.56   4.19 1.02      395      470
##      m      2.26   2.25 0.21 0.20   1.95   2.63 1.00     1919     1555
##      s      0.94   0.91 0.16 0.15   0.72   1.25 1.01      441      622
##      y1     2.26   2.13 0.95 0.87   0.96   3.96 1.00     1786     1868

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

log normal dist.

logalithm of variable follows normal distribution
log y~N(m,s), y>0
E[y]=exp(m+s^2)
V[y]=exp(2*m+s^2)*(exp(s^2)-1)

ex1-7.stan

data{
  int N;
  vector[N]<lower=0> y;
}
parameters{
  real m0;
  real<lower=0> s0;
}
model{
  y~lognormal(m0,s0);
}
n=30
y=rlnorm(n,0.5,1)
summary(y)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.39    0.95    2.05    3.26    4.67   14.38
par(mfrow=c(1,2))
hist(y,main='y')
density(y) |> plot(main='y')

data=list(N=n,y=y)

mdl=cmdstan_model('./ex1-7.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -77.4717 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        7       -41.985   8.65034e-05   0.000161723      0.9434      0.9434       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__   -41.98
##      m0       0.73
##      s0       0.98
##      m        5.45
##      s        4.28
##      y1       7.25
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable   mean median   sd  mad     q5    q95 rhat ess_bulk ess_tail
##      lp__ -43.03 -42.74 1.00 0.75 -45.07 -42.05 1.00      921     1444
##      m0     0.74   0.75 0.19 0.19   0.42   1.04 1.00      966      971
##      s0     1.04   1.03 0.15 0.14   0.83   1.30 1.00     1917     1437
##      m      6.86   6.07 3.42 2.11   3.74  12.29 1.00     1250     1166
##      s      5.68   4.86 3.35 2.01   2.69  11.04 1.00     1318     1175
##      y1     3.47   1.98 4.90 1.77   0.38  12.15 1.00     2031     1958

y1=mcmc$draws('y1')
par(mfrow=c(1,2))
hist(y1,main='y1')
density(y1) |> plot(main='y1') 

beta dist.

use as prior of binomial dist parameter p
event with probabilty p0 occur a-1 times, do not occur b-1 times in a+b-2 trials
p~Be(a,b), p[0,1], p=p0^(a-1)*(1-p0)^(b-1)
E[p]=a/(a+b), V[p]=ab/(a+b)^2/(a+b+1)

ex1-8.stan

data{
  int N;
  vector[N] p;
}
parameters{
  real<lower=0> a;
  real<lower=0> b;
}
model{
  p~beta(a,b);
}
n=20
p=rbeta(n,1,2)
summary(p)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.007   0.114   0.264   0.333   0.446   0.961
par(mfrow=c(1,2))
hist(p,main='p')
density(p) |> plot(main='p')

data=list(N=n,p=p)

mdl=cmdstan_model('./ex1-8.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = -184.562 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##       11       3.13888   7.26822e-05   7.60616e-07           1           1       12    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##      lp__     3.14
##      a        0.76
##      b        1.42
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable mean median   sd  mad   q5  q95 rhat ess_bulk ess_tail
##      lp__ 2.40   2.69 0.95 0.68 0.43 3.30 1.00      906     1305
##      a    0.86   0.84 0.22 0.21 0.53 1.25 1.01      701     1154
##      b    1.66   1.61 0.47 0.46 0.97 2.50 1.01      632      877

Dirichlet dist.

use as prior of categorical dist parameter p
p~dir(p0), p[0,1]

ex1-9.stan

data {
  int N;
  int K;
  matrix[N, K] p;
}
parameters {
  simplex[K] p0;
}
model {
  for(i in 1:N){
     p[i]~dirichlet(p0);
  }
}
library(gtools)
n=20
p=rdirichlet(n,c(0.5,0.3,0.2))
summary(p)
##        V1              V2              V3       
##  Min.   :0.001   Min.   :0.000   Min.   :0.000  
##  1st Qu.:0.257   1st Qu.:0.028   1st Qu.:0.000  
##  Median :0.827   Median :0.066   Median :0.013  
##  Mean   :0.622   Mean   :0.258   Mean   :0.120  
##  3rd Qu.:0.944   3rd Qu.:0.456   3rd Qu.:0.105  
##  Max.   :0.999   Max.   :0.945   Max.   :0.996
boxplot(p)

data=list(N=n,K=ncol(p),p=p)

mdl=cmdstan_model('./ex1-9.stan')

mle=mdl$optimize(data=data)
## Initial log joint probability = 51.4479 
##     Iter      log prob        ||dx||      ||grad||       alpha      alpha0  # evals  Notes  
##        6       94.4034   0.000122507   0.000508777           1           1       10    
## Optimization terminated normally:  
##   Convergence detected: relative gradient magnitude is below tolerance 
## Finished in  0.1 seconds.
mle
##  variable estimate
##     lp__     94.40
##     p0[1]     0.56
##     p0[2]     0.29
##     p0[3]     0.15
mcmc=goMCMC(mdl,data)
## Running MCMC with 4 parallel chains...
## 
## Chain 1 WARNING: There aren't enough warmup iterations to fit the 
## Chain 1          three stages of adaptation as currently configured. 
## Chain 1          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 1          the given number of warmup iterations: 
## Chain 1            init_buffer = 15 
## Chain 1            adapt_window = 75 
## Chain 1            term_buffer = 10 
## Chain 1 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 1 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 1 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 2 WARNING: There aren't enough warmup iterations to fit the 
## Chain 2          three stages of adaptation as currently configured. 
## Chain 2          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 2          the given number of warmup iterations: 
## Chain 2            init_buffer = 15 
## Chain 2            adapt_window = 75 
## Chain 2            term_buffer = 10 
## Chain 2 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 2 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 2 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 3 WARNING: There aren't enough warmup iterations to fit the 
## Chain 3          three stages of adaptation as currently configured. 
## Chain 3          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 3          the given number of warmup iterations: 
## Chain 3            init_buffer = 15 
## Chain 3            adapt_window = 75 
## Chain 3            term_buffer = 10 
## Chain 3 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 3 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 3 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 4 WARNING: There aren't enough warmup iterations to fit the 
## Chain 4          three stages of adaptation as currently configured. 
## Chain 4          Reducing each adaptation stage to 15%/75%/10% of 
## Chain 4          the given number of warmup iterations: 
## Chain 4            init_buffer = 15 
## Chain 4            adapt_window = 75 
## Chain 4            term_buffer = 10 
## Chain 4 Iteration:   1 / 600 [  0%]  (Warmup) 
## Chain 4 Iteration: 101 / 600 [ 16%]  (Sampling) 
## Chain 4 Iteration: 600 / 600 [100%]  (Sampling) 
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
seeMCMC(mcmc)
##  variable  mean median   sd  mad    q5   q95 rhat ess_bulk ess_tail
##     lp__  89.69  90.01 1.06 0.72 87.68 90.64 1.01      928     1005
##     p0[1]  0.55   0.55 0.06 0.05  0.45  0.64 1.00     1933     1392
##     p0[2]  0.30   0.30 0.05 0.05  0.21  0.39 1.00     1676     1289
##     p0[3]  0.15   0.15 0.03 0.03  0.11  0.21 1.00     1873     1494